#include <src/pcps.h>
#include <src/input.h>
#include <src/correlator.h>
#include <src/wavefunction.h>
#include <src/hamiltonian.h>
#include <src/cluster.h>
#include <src/random.h>
#include <src/blas_lapack.h>
#include <src/hf.h>

#include <ctime>
#include <stack>
#include <algorithm>

//-------------------------------------------------------------------------------
// correlators_from_site_sets
//                                                  
//                          create correlators from sets of sites
//-------------------------------------------------------------------------------

template <class S>
inline void correlators_from_site_sets(const pcps::Input & userinp,
                                       const std::set< std::set<int> > & site_sets,
                                       std::vector< pcps::Correlator<S> > & correlators) {

  // create a correlator from each set of sites
  for (std::set< std::set<int> >::const_iterator sites = site_sets.begin(); sites != site_sets.end(); sites++) {

    // skip empty sets
    if (sites->size() == 0) {

    // create fermion correlators
    } else if ( userinp.particle_type() == pcps::keyword::Part_Fermion ) {

      std::set<int> alpha_and_beta_sites;
      for (std::set<int>::const_iterator site = sites->begin(); site != sites->end(); site++) {
        alpha_and_beta_sites.insert( 2 * (*site)    ); // alpha site
        alpha_and_beta_sites.insert( 2 * (*site) + 1); //  beta site
      }
      correlators.push_back( pcps::Correlator<S>(alpha_and_beta_sites, userinp.configs_per_site(), pcps::unity<S>()) );

    // create correlators for other particle types
    } else correlators.push_back( pcps::Correlator<S>(*sites, userinp.configs_per_site(), pcps::unity<S>()) );

  }

}

//-------------------------------------------------------------------------------
// sign_rule_guess -- appends one-site correlators to the inputted vector
//                    of correlators that establish Marshall's sign rule
//-------------------------------------------------------------------------------

template <class S>
inline void sign_rule_guess(const pcps::Input &userinp,
                            std::vector< pcps::Correlator<S> > &correlators) {

  // check that the lattice is cartesian
  if ( userinp.lattice_type() != pcps::keyword::Lattice_Cartesian )
    throw pcps::Exception( (boost::format("sign rule guess only implemented for lattice_type \"%s\"") % pcps::keyword::Lattice_Cartesian).str() );

  // check that we are dealing with spin-half spins
  if (userinp.particle_type() != pcps::keyword::Part_Spin || userinp.configs_per_site() != 2)
    throw pcps::Exception( (boost::format("sign rule guess requires spin-half spins")).str() );

  // create a function object that will recursively find the signed sites
  struct Local_Func_Object {
    static void get_signed_sites(const int lx,
                                 const int ly,
                                 const int lz,
                                 std::set<int> &sites,
                                 const int ix = 0,
                                 const int iy = 0,
                                 const int iz = 0,
                                 const bool last_was_A = false) {

      // end of recursion
      if (ix == lx || iy == ly || iz == lz)
        return;

      // check whether this site is an A site
      const bool this_is_A = !last_was_A;

      // if this site is a B site, add it to the list of signed sites
      if (!this_is_A)
        sites.insert( lx * ly * iz + lx * iy + ix );

      // recursion for next sites
      get_signed_sites(lx, ly, lz, sites, ix+1, iy, iz, this_is_A);
      get_signed_sites(lx, ly, lz, sites, ix, iy+1, iz, this_is_A);
      get_signed_sites(lx, ly, lz, sites, ix, iy, iz+1, this_is_A);

    }
  };

  // get the lattice dimensions
  const int lx = userinp.lattice_length_x();
  const int ly = userinp.lattice_length_y();
  const int lz = userinp.lattice_length_z();

  // get the set of sites which should create signs when their spin is down
  std::set<int> signed_sites;
  Local_Func_Object::get_signed_sites(lx, ly, lz, signed_sites);

  // determine whether or not to fix the sign rule correlators
  const bool fixed = userinp.fix_classical();

  // for each signed site, create the necessary signs in the wavefunction by adding a one-site correlator
  for (std::set<int>::const_iterator site = signed_sites.begin(); site != signed_sites.end(); site++) {
    pcps::Correlator<S> corr(*site, userinp.configs_per_site(), pcps::unity<S>(), fixed);
    corr[1] *= -1.0;
    correlators.push_back(corr);
  }

}

//-------------------------------------------------------------------------------
// pcps::get_classical_tripartite_guess -- construct an initial guess for a
//                                         tripartite lattice in which each
//                                         sublattice's spin's expectation value
//                                         is rotated at an angle of 2*pi/3
//                                         (120 degrees) from the others
//-------------------------------------------------------------------------------

template <class S>
void pcps::get_classical_tripartite_guess(const pcps::Input &userinp,
                                          const std::set<int> & a_site_set,
                                          const std::set<int> & b_site_set,
                                          const std::set<int> & c_site_set,
                                          std::vector< pcps::Correlator<S> > &correlators) {

  // check that we are using complex numbers
  if (typeid(S) != typeid(std::complex<double>))
    throw pcps::Exception( (boost::format("pcps::get_classical_tripartite_guess requires complex numbers")).str() );

  // check that we are dealing with spins
  if (userinp.particle_type() != pcps::keyword::Part_Spin)
    throw pcps::Exception( (boost::format("pcps::get_classical_tripartite_guess requires particles to be spins")).str() );

  const double pi = 3.14159265358979;

  const double x = 0.0 * pi;

  // determine whether or not to fix the classical correlators
  const bool fixed = userinp.fix_classical();

  // a sublattice
  {
    const double theta = 1.0 * pi / 2.0;
    const double phi = 0.0 * pi / 3.0;
    for (std::set<int>::const_iterator site = a_site_set.begin(); site != a_site_set.end(); site++)
      correlators.push_back( pcps::spin_coherent_correlator<S>(*site, userinp.configs_per_site(), theta, phi, fixed) );
  }

  // b sublattice
  {
    const double theta = 1.0 * pi / 2.0 + x;
    const double phi = 2.0 * pi / 3.0;
    for (std::set<int>::const_iterator site = b_site_set.begin(); site != b_site_set.end(); site++)
      correlators.push_back( pcps::spin_coherent_correlator<S>(*site, userinp.configs_per_site(), theta, phi, fixed) );
  }

  // c sublattice
  {
    const double theta = 1.0 * pi / 2.0 - x;
    const double phi = 4.0 * pi / 3.0;
    for (std::set<int>::const_iterator site = c_site_set.begin(); site != c_site_set.end(); site++)
      correlators.push_back( pcps::spin_coherent_correlator<S>(*site, userinp.configs_per_site(), theta, phi, fixed) );
  }

}

//-------------------------------------------------------------------------------
// pcps::get_classical_icosidodecahedron_guess -- construct an initial guess
//                                                using single site correlators
//                                                to express the spin coherent
//                                                state representation of the
//                                                classical heisenberg solution
//-------------------------------------------------------------------------------

template <class S>
void pcps::get_classical_icosidodecahedron_guess(const pcps::Input &userinp, std::vector< pcps::Correlator<S> > &correlators) {

  // check that we are using complex numbers
  if (typeid(S) != typeid(std::complex<double>))
    throw pcps::Exception( (boost::format("pcps::get_classical_icosidodecahedron_guess requires complex numbers")).str() );

  // check that we are dealing with spins
  if (userinp.particle_type() != pcps::keyword::Part_Spin)
    throw pcps::Exception( (boost::format("pcps::get_classical_icosidodecahedron_guess requires particles to be spins")).str() );

  // check that the lattice is the icosidodecahedron
  if ( userinp.lattice_type() != pcps::keyword::Lattice_Icosidodecahedron )
    throw pcps::Exception( (boost::format("pcps::get_classical_icosidodecahedron_guess requires the icosidodecahedron lattice")).str() );

  int a_site_array [10] = { 27, 21, 23,  1,  2, 17,  7, 11, 15,  6 };
  int b_site_array [10] = { 24, 28, 22, 26,  3,  5, 14, 16, 10, 12 };
  int c_site_array [10] = { 30, 25, 29, 19, 20,  4,  8, 18, 13,  9 };

  // count from zero
  for (int i = 0; i < 10; i++) {
    a_site_array [i] -= 1;
    b_site_array [i] -= 1;
    c_site_array [i] -= 1;
  }

  // create sets
  std::set<int> a_site_set(&a_site_array[0], &a_site_array[0] + 10);
  std::set<int> b_site_set(&b_site_array[0], &b_site_array[0] + 10);
  std::set<int> c_site_set(&c_site_array[0], &c_site_array[0] + 10);
  assert(a_site_set.size() == 10);
  assert(b_site_set.size() == 10);
  assert(c_site_set.size() == 10);

  // create the correlators
  pcps::get_classical_tripartite_guess(userinp, a_site_set, b_site_set, c_site_set, correlators);

}

//-------------------------------------------------------------------------------
// pcps::get_classical_triangular_guess -- construct an initial guess using
//                                         single site correlators to express the
//                                         spin coherent state representation of
//                                         the classical heisenberg solution of
//                                         the triangular lattice
//-------------------------------------------------------------------------------

template <class S>
void pcps::get_classical_triangular_guess(const pcps::Input &userinp, std::vector< pcps::Correlator<S> > &correlators) {

  // check that we are using complex numbers
  if (typeid(S) != typeid(std::complex<double>))
    throw pcps::Exception( (boost::format("pcps::get_classical_triangular_guess requires complex numbers")).str() );

  // check that we are dealing with spins
  if (userinp.particle_type() != pcps::keyword::Part_Spin)
    throw pcps::Exception( (boost::format("pcps::get_classical_triangular_guess requires particles to be spins")).str() );

  // check that the lattice is triangular
  if ( userinp.lattice_type() != pcps::keyword::Lattice_Triangular )
    throw pcps::Exception( (boost::format("pcps::get_classical_triangular_guess requires the icosidodecahedron lattice")).str() );

  // get the lattice parameters
  const int length_x = userinp.lattice_length_x();
  const int length_y = userinp.lattice_length_y();
  const int length_z = userinp.lattice_length_z();

  // check that the lattice is 2 dimensional
  if (length_z != 1)
    throw pcps::Exception("pcps::get_classical_triangular_guess requires lattice_length_z to equal 1");

  // check that no lengths are negative
  if (length_x <= 0 || length_y <= 0 || length_z <= 0)
    throw pcps::Exception("pcps::get_classical_triangular_guess requires lattice lengths to be positive integers");

  // create sets of sites for each sublattice
  std::set<int> a_site_set;
  std::set<int> b_site_set;
  std::set<int> c_site_set;

  // assign each lattice site to a sublattice
  for (int i = 0; i < length_x; i++) 
  for (int j = 0; j < length_y; j++) 
  for (int k = 0; k < length_z; k++) {
    if ( (i+j) % 3 == 0 ) a_site_set.insert(length_x * length_y * k + length_x * j + i);
    if ( (i+j) % 3 == 1 ) b_site_set.insert(length_x * length_y * k + length_x * j + i);
    if ( (i+j) % 3 == 2 ) c_site_set.insert(length_x * length_y * k + length_x * j + i);
  }

  // create the correlators
  pcps::get_classical_tripartite_guess(userinp, a_site_set, b_site_set, c_site_set, correlators);

}

//-------------------------------------------------------------------------------
// pcps::Wavefunction::get_text_correlators
//
//                         Reads correlators from a file and appends them to the
//                         wavefunction's vector of correlators.
//-------------------------------------------------------------------------------

template <class S>
void pcps::Wavefunction<S>::get_text_correlators(const pcps::Input & userinp) {

  // allocate a buffer for storing file lines
  const int buffer_size = 4096;
  boost::shared_array<char> buffer = pcps::allocate_shared_array<char>(buffer_size+1);

  // get the name of the file holding the correlator data
  const std::string filename = userinp.workdir() + "pcps_correlators.txt";

  // open the file for reading
  std::ifstream ifs(filename.c_str());

  // check that the file opened successfully
  if (!ifs.is_open())
    throw pcps::Exception((boost::format("failed to open file \"%s\" in pcps::Wavefunction::get_text_correlators") % filename).str());

  // process the file line by line
  while (!ifs.eof()) {

    // get the next line of the file
    ifs.getline(buffer.get(), buffer_size);
    std::string line(buffer.get());

    // remove comments
    pcps::strip_comments(line, '%');

    // split the line into words
    std::vector<std::string> word_vec = pcps::split_string(line, " ");

    // do nothing if the line is empty
    if (word_vec.size() == 0) {}

    // read a correlator
    else if (word_vec.at(0) == "correlator") {

      std::vector<int> sites, ti_sites;

      bool have_sites = false;
      bool fixed = false;
      bool ti = false;

      std::vector< std::vector<std::string> > value_lines;

      while (true) {

        // make sure we are not at the end of the file
        if (ifs.eof())
          throw pcps::Exception("error in pcps::Wavefunction::get_text_correlators.  Encountered end of file while reading a correlator.");

        // get the next line of the file
        ifs.getline(buffer.get(), buffer_size);
        std::string line(buffer.get());

        // remove comments
        pcps::strip_comments(line, '%');

        // split the line into words
        std::vector<std::string> word_vec = pcps::split_string(line, " ");

        if (word_vec.size() == 0) {}

        else if (word_vec.at(0) == "sites") {
          if (have_sites) throw pcps::Exception("error in pcps::Wavefunction::get_text_correlators.  Each correlator may have only one \"sites\" line.");
          have_sites = true;
          for (int i = 1; i < word_vec.size(); i++)
            sites.push_back( std::atoi(word_vec.at(i).c_str()) );
        }

        else if (word_vec.at(0) == "ti_sites") {
          if (ti) throw pcps::Exception("error in pcps::Wavefunction::get_text_correlators.  Each correlator may have only one \"ti_sites\" line.");
          ti = true;
          for (int i = 1; i < word_vec.size(); i++)
            ti_sites.push_back( std::atoi(word_vec.at(i).c_str()) );
        }

        else if (word_vec.at(0) == "fixed")
          fixed = true;

        else if (word_vec.at(0) == "end")
          break;

        else
          value_lines.push_back(word_vec);

      }

      if (ti) {

        if (fixed) throw pcps::Exception("error in pcps::Wavefunction::get_text_correlators.  ti correlator cannot be specified as fixed.  Whether it is fixed is determined by the base correlator it associates with.");

        if (value_lines.size() > 0) throw pcps::Exception("error in pcps::Wavefunction::get_text_correlators.  ti correlator cannot have value lines.");

        // create the translationally invariant correlator
        pcps::Correlator<S> ti_corr(_correlators.at(this->get_ti_parent_index(ti_sites)), sites);
        _correlators.push_back(ti_corr);

      } else {

        // create a regular correlator
        _correlators.push_back( pcps::Correlator<S>(sites, userinp.configs_per_site(), pcps::unity<S>(), fixed) );

      }

      // process value lines, which specify specific values for some of the correlator's elements
      for ( std::vector< std::vector<std::string> >::const_iterator words = value_lines.begin(); words != value_lines.end(); words++) {

        const int nval = ( typeid(S) == typeid(std::complex<double>) ? 2 : 1 );

        if (words->size() != _correlators.rbegin()->site_vec().size() + nval)
          throw pcps::Exception("error in pcps::Wavefunction::get_text_correlators. Bad value line.");

        int cmpd = 0;
        for (int j = 0; j < _correlators.rbegin()->site_vec().size(); j++) {
          int occ = std::atoi(words->at(j).c_str());
          cmpd = cmpd * userinp.configs_per_site() + occ;
        }
        (*_correlators.rbegin())[cmpd] = pcps::unity<S>() * std::atof(words->at(_correlators.rbegin()->site_vec().size()).c_str());
        if ( typeid(S) == typeid(std::complex<double>) )
          (*_correlators.rbegin())[cmpd] += pcps::imaginary_unity<S>() * std::atof(words->at(_correlators.rbegin()->site_vec().size()+1).c_str());

      }

    }

    // unknown keyword
    else
      throw pcps::Exception( ( boost::format("unknown keyword \"%s\" when reading file \"%s\"") % word_vec.at(0) % filename ).str() );

  }

  // close the file
  ifs.close();

}

//-------------------------------------------------------------------------------
// pcps::Wavefunction::get_binary_correlators -- read correlators from a
//                                               binary file
//-------------------------------------------------------------------------------

template <class S>
void pcps::Wavefunction<S>::get_binary_correlators(const pcps::Input & userinp) {

  // get the name of the binary file
  const std::string filename = userinp.workdir() + "pcps_wfn.bin";

  // open the binary file for reading
  std::ifstream ifs(filename.c_str(), std::ios::in | std::ios::binary);

  // check that the file is open
  if (!ifs.is_open())
    throw pcps::Exception( (boost::format("failed to open file \"%s\" in pcps::Wavefunction::get_binary_correlators") % filename).str() );

  // read the number of correlators
  int n;
  ifs.read((char *)&n, sizeof(int));

  // read the correlators
  for (int i = 0; i < n; i++) {
    pcps::Correlator<S> new_corr;
    new_corr.read_binary(ifs);
    _correlators.push_back(new_corr);
  }

  // close the file
  ifs.close();

}

//-------------------------------------------------------------------------------
// pcps::Wavefunction::get_classical_correlators
//                                                  
//                       create correlators corresponding to a classical state
//-------------------------------------------------------------------------------

template <class S>
void pcps::Wavefunction<S>::get_classical_correlators(const pcps::Input & userinp) {

  // for spins on a cartesian lattice, use Marshall's sign rule as a guess
  if      (    userinp.particle_type() == pcps::keyword::Part_Spin
            && userinp.lattice_type() == pcps::keyword::Lattice_Cartesian )
    sign_rule_guess(userinp, _correlators);

  // for spins on the icosidodecahedron lattice, use the solution to the classical heisenberg model as a guess
  else if (    userinp.particle_type() == pcps::keyword::Part_Spin
            && userinp.lattice_type() == pcps::keyword::Lattice_Icosidodecahedron )
    pcps::get_classical_icosidodecahedron_guess(userinp, _correlators);

  // for spins on the triangular lattice, use the solution to the classical heisenberg model as a guess
  else if (    userinp.particle_type() == pcps::keyword::Part_Spin
            && userinp.lattice_type() == pcps::keyword::Lattice_Triangular )
    pcps::get_classical_triangular_guess(userinp, _correlators);

  // throw an exception if there is no classical guess implemented
  else
    throw pcps::Exception( (boost::format("no classical guess implemented for the specified particle and lattice types")).str() );

}

//-------------------------------------------------------------------------------
// pcps::Wavefunction::get_lattice_correlators
//                                                  
//                        create correlators based on the lattice
//-------------------------------------------------------------------------------

template <class S>
void pcps::Wavefunction<S>::get_lattice_correlators(const pcps::Input & userinp) {

  // initialize the sets of sites for the correlators
  std::set< std::set<int> > corr_site_sets;

  // create the correlators' site sets on the cartesian lattice
  if ( userinp.lattice_type() == pcps::keyword::Lattice_Cartesian ) {

    // get the lattice dimensions
    const int lx = userinp.lattice_length_x();
    const int ly = userinp.lattice_length_y();
    const int lz = userinp.lattice_length_z();

    // get correlator size and boundary type
    const int corr_size = userinp.correlator_size();
    const bool periodic = ( userinp.boundary() == pcps::keyword::Boundary_Periodic );

    // loop over all sites on the lattice
    for (int i = 0; i < lx; i++) 
    for (int j = 0; j < ly; j++) 
    for (int k = 0; k < lz; k++) {

      // create line correlators anchored at this site
      if (userinp.correlator_type() == pcps::keyword::Correlator_Line) {

        // if there is room to make a line correlator couting upwards from the site along the x axis, do so
        if (periodic || i - 1 + corr_size < lx) {
          std::set<int> temp_set;
          for (int p = 0; p < corr_size; p++)
            temp_set.insert( lx * ly * k + lx * j + (i+p) % lx );
          if (temp_set.size() == corr_size)
            corr_site_sets.insert(temp_set);
        }

        // if there is room to make a line correlator couting upwards from the site along the y axis, do so
        if (periodic || j - 1 + corr_size < ly) {
          std::set<int> temp_set;
          for (int p = 0; p < corr_size; p++)
            temp_set.insert( lx * ly * k + lx * ( (j+p) % ly ) + i );
          if (temp_set.size() == corr_size)
            corr_site_sets.insert(temp_set);
        }

        // if there is room to make a line correlator couting upwards from the site along the z axis, do so
        if (periodic || k - 1 + corr_size < lz) {
          std::set<int> temp_set;
          for (int p = 0; p < corr_size; p++)
            temp_set.insert( lx * ly * ( (k+p) % lz ) + lx * j + i );
          if (temp_set.size() == corr_size)
            corr_site_sets.insert(temp_set);
        }

      // create square correlators anchored at this site
      } else if (userinp.correlator_type() == pcps::keyword::Correlator_Square) {

        // if there is room to make a square correlator couting upwards from the site along the x and y axes, do so
        if ( periodic || ( i - 1 + corr_size < lx && j - 1 + corr_size < ly ) ) {
          std::set<int> temp_set;
          for (int p = 0; p < corr_size; p++)
          for (int q = 0; q < corr_size; q++)
            temp_set.insert( lx * ly * k + lx * ( (j+q) % ly ) + (i+p) % lx );
          if ( temp_set.size() == corr_size * corr_size )
            corr_site_sets.insert(temp_set);
        }

        // if there is room to make a square correlator couting upwards from the site along the x and z axes, do so
        if ( periodic || ( i - 1 + corr_size < lx && k - 1 + corr_size < lz ) ) {
          std::set<int> temp_set;
          for (int p = 0; p < corr_size; p++)
          for (int q = 0; q < corr_size; q++)
            temp_set.insert( lx * ly * ( (k+q) % lz ) + lx * j + (i+p) % lx );
          if ( temp_set.size() == corr_size * corr_size )
            corr_site_sets.insert(temp_set);
        }

        // if there is room to make a square correlator couting upwards from the site along the y and z axes, do so
        if ( periodic || ( j - 1 + corr_size < ly && k - 1 + corr_size < lz ) ) {
          std::set<int> temp_set;
          for (int p = 0; p < corr_size; p++)
          for (int q = 0; q < corr_size; q++)
            temp_set.insert( lx * ly * ( (k+q) % lz ) + lx * ( (j+p) % ly ) + i );
          if ( temp_set.size() == corr_size * corr_size )
            corr_site_sets.insert(temp_set);
        }

      // unknown correlator type
      } else throw pcps::Exception( (boost::format("unknown correlator type for cartesian lattice: \"%s\"") % userinp.correlator_type()).str() );

    }

  // create the correlators' site sets on the triangular lattice
  } else if ( userinp.lattice_type() == pcps::keyword::Lattice_Triangular ) {

    // get the lattice dimensions
    const int lx = userinp.lattice_length_x();
    const int ly = userinp.lattice_length_y();
    const int lz = userinp.lattice_length_z();

    // check that lattice is not 3-dimensional
    if (lz != 1) throw pcps::Exception( (boost::format("lattice_length_z must be 1 for the triangular lattice")).str() );

    // get correlator size and boundary type
    const int corr_size = userinp.correlator_size();
    const bool periodic = ( userinp.boundary() == pcps::keyword::Boundary_Periodic );

    // loop over all sites on the lattice
    for (int i = 0; i < lx; i++)
    for (int j = 0; j < ly; j++)
    for (int k = 0; k < lz; k++) {

      // create line correlators anchored at this site
      if (userinp.correlator_type() == pcps::keyword::Correlator_Line) {

        // if there is room to make a line correlator along the a direction, do so
        if (periodic || i - 1 + corr_size < lx) {
          std::set<int> temp_set;
          for (int p = 0; p < corr_size; p++)
            temp_set.insert( lx * ly * k + lx * j + (i+p) % lx );
          if (temp_set.size() == corr_size)
            corr_site_sets.insert(temp_set);
        }

        // if there is room to make a line correlator along the b direction, do so
        if (periodic || j - 1 + corr_size < ly) {
          std::set<int> temp_set;
          for (int p = 0; p < corr_size; p++)
            temp_set.insert( lx * ly * k + lx * ( (j+p) % ly ) +  i );
          if (temp_set.size() == corr_size)
            corr_site_sets.insert(temp_set);
        }

        // if there is room to make a line correlator along the c direction, do so
        if (periodic || ( i - 1 + corr_size < lx && j - 1 + corr_size < ly ) ) {
          std::set<int> temp_set;
          for (int p = 0; p < corr_size; p++)
            temp_set.insert( lx * ly * k + lx * ( (j+p) % ly ) + (i+p) % lx );
          if (temp_set.size() == corr_size)
            corr_site_sets.insert(temp_set);
        }

      // create triangle correlators anchored at this site
      } else if (userinp.correlator_type() == pcps::keyword::Correlator_Triangle) {

        // if there is room to make a triangle correlator in the (+0,+1), (+1,+1) direction, do so
        if ( periodic || ( i - 1 + corr_size < lx && j - 1 + corr_size < ly ) ) {
          std::set<int> temp_set;
          for (int q = 0; q < corr_size; q++)
          for (int p = 0; p <= q; p++)
            temp_set.insert( lx * ly * k + lx * ( (j+q) % ly ) + (i+p) % lx );
          if ( temp_set.size() == corr_size * (corr_size+1) / 2 )
            corr_site_sets.insert(temp_set);
        }

        // if there is room to make a triangle correlator in the (+0,-1), (-1,-1) direction, do so
        if ( periodic || ( i + 1 - corr_size >= 0 && j + 1 - corr_size >= 0 ) ) {
          std::set<int> temp_set;
          for (int q = 0; q < corr_size; q++)
          for (int p = 0; p <= q; p++)
            temp_set.insert( lx * ly * k + lx * ( (j-q+ly) % ly ) + (i-p+lx) % lx );
          if ( temp_set.size() == corr_size * (corr_size+1) / 2 )
            corr_site_sets.insert(temp_set);
        }

      // create square (rhombus) correlators anchored at this site
      } else if (userinp.correlator_type() == pcps::keyword::Correlator_Square) {

        // if there is room to make a rhombus correlator in the (+1,-1) direction, do so
        if ( periodic || ( i - 1 + corr_size < lx && j + 1 - corr_size >= 0 ) ) {
          std::set<int> temp_set;
          for (int p = 0; p < corr_size; p++)
          for (int q = 0; q < corr_size; q++)
            temp_set.insert( lx * ly * k + lx * ( (j-q+ly) % ly ) + (i+p) % lx );
          if ( temp_set.size() == corr_size * corr_size )
            corr_site_sets.insert(temp_set);
        }

        // if there is room to make a rhombus correlator in the (+1,+2) direction, do so
        if ( periodic || ( i - 1 + corr_size < lx && j + 2*(corr_size-1) < ly ) ) {
          std::set<int> temp_set;
          for (int p = 0; p < corr_size; p++)
          for (int q = 0; q < corr_size; q++)
            temp_set.insert( lx * ly * k + lx * ( (j+q+p) % ly ) + (i+p) % lx );
          if ( temp_set.size() == corr_size * corr_size )
            corr_site_sets.insert(temp_set);
        }

        // if there is room to make a rhombus correlator in the (-2,-1) direction, do so
        if ( periodic || ( i - 2*(corr_size-1) >= 0 && j - (corr_size-1) >= 0 ) ) {
          std::set<int> temp_set;
          for (int p = 0; p < corr_size; p++)
          for (int q = 0; q < corr_size; q++)
            temp_set.insert( lx * ly * k + lx * ( (j-q+ly) % ly ) + (i-p-q+lx+lx) % lx );
          if ( temp_set.size() == corr_size * corr_size )
            corr_site_sets.insert(temp_set);
        }

      // unknown correlator type
      } else throw pcps::Exception( (boost::format("unknown correlator type for triangular lattice: \"%s\"") % userinp.correlator_type()).str() );

    }

  // throw an error if the lattice type is not supported
  } else throw pcps::Exception((boost::format("pcps::Wavefunction::get_lattice_correlators not defined for \"%s\" lattice") % userinp.lattice_type()).str());

  // create correlators from the sets of sites
  correlators_from_site_sets(userinp, corr_site_sets, _correlators);

}

//-------------------------------------------------------------------------------
// pcps::Wavefunction::get_all_pair_correlators
//                                                  
//                        create a correlator for each pair of sites
//-------------------------------------------------------------------------------

template <class S>
void pcps::Wavefunction<S>::get_all_pair_correlators(const pcps::Input &userinp, std::vector<pcps::Correlator<S> > & correlators) {

  // determine the number of sites to pair
  const int n = ( userinp.particle_type() == pcps::keyword::Part_Fermion ? userinp.nsites() / 2 : userinp.nsites() );

  // get the site pairs
  std::set< std::set<int> > site_pairs;
  for (int i = 0; i < n; i++)
  for (int j = i+1; j < n; j++) {
    std::set<int> s;
    s.insert(i);
    s.insert(j);
    site_pairs.insert(s);
  }

  // create correlators from the site pairs
  correlators_from_site_sets(userinp, site_pairs, correlators);

}

//-------------------------------------------------------------------------------
// pcps::Wavefunction::init_correlators -- initialize the correlators
//-------------------------------------------------------------------------------

template <class S>
void pcps::Wavefunction<S>::init_correlators(const pcps::Input & userinp) {

  // get MPI info
  const MPI::Comm & comm = MPI::COMM_WORLD;
  const int nproc = comm.Get_size();
  const int myrank = comm.Get_rank();

  // clear the vectors of correlators
  _correlators.clear();
  _pert_corrs.clear();

  // prepare correlators on root node
  if (myrank == 0) {

    // check that all correlator sources are recognized
    for (std::set<std::string>::const_iterator source = userinp.corr_sources().begin(); source != userinp.corr_sources().end(); source++) {
      if (    *source != pcps::keyword::Corr_Source_None
           && *source != pcps::keyword::Corr_Source_Binary
           && *source != pcps::keyword::Corr_Source_Text
           && *source != pcps::keyword::Corr_Source_Classical
           && *source != pcps::keyword::Corr_Source_Lattice
           && *source != pcps::keyword::Corr_Source_All_Pair )
        throw pcps::Exception( (boost::format("unknown correlator source \"%s\" in pcps::Wavefunction::init_correlators") % *source).str() );
    }

    // for each specified correlator source, get the appropriate correlators
    if ( userinp.corr_sources().count( pcps::keyword::Corr_Source_Binary    ) > 0 ) this->get_binary_correlators(userinp);
    if ( userinp.corr_sources().count( pcps::keyword::Corr_Source_Lattice   ) > 0 ) this->get_lattice_correlators(userinp);
    if ( userinp.corr_sources().count( pcps::keyword::Corr_Source_All_Pair  ) > 0 ) this->get_all_pair_correlators(userinp, _correlators);
    if ( userinp.corr_sources().count( pcps::keyword::Corr_Source_Classical ) > 0 ) this->get_classical_correlators(userinp);
    if ( userinp.corr_sources().count( pcps::keyword::Corr_Source_Text      ) > 0 ) this->get_text_correlators(userinp);

    // absorb any non-fixed correlators contained by other correlators
    pcps::absorb_correlators(userinp, _correlators);

    // if requested, set non-fixed and non-ti correlator elements to have unit norm and random phase
    if (userinp.random_correlators()) {
      for (typename std::vector<pcps::Correlator<S> >::iterator corr = _correlators.begin(); corr != _correlators.end(); corr++) {
        if (corr->fixed() || corr->ti()) continue;
        for (int i = 0; i < corr->nelements(); i++) {
          do {
            (*corr)[i] =             pcps::unity<S>() * ( pcps::random_number<double>() - 0.5 )
                         + pcps::imaginary_unity<S>() * ( pcps::random_number<double>() - 0.5 );
          } while ( std::abs((*corr)[i]) < 1.0e-4 );
          (*corr)[i] /= std::abs((*corr)[i]);
        }
      }
    }

    // create perturbative pair-correlators (omit those contained by non-perturbative correlators)
    if (userinp.do_perturb()) {
      this->get_all_pair_correlators(userinp, _pert_corrs);
      for (int i = _pert_corrs.size()-1; i >= 0; i--) {
        bool contained = false;
        for (typename std::vector< pcps::Correlator<S> >::const_iterator corr = _correlators.begin(); corr != _correlators.end(); corr++) {
          int n_contained = 0;
          for (std::vector<int>::const_iterator site = _pert_corrs.at(i).site_vec().begin(); site != _pert_corrs.at(i).site_vec().end(); site++)
            n_contained += std::count(corr->site_vec().begin(), corr->site_vec().end(), *site);
          if ( n_contained == _pert_corrs.at(i).site_vec().size() )
            contained = true;
        }
        if (contained)
          _pert_corrs.erase( _pert_corrs.begin() + i );
      }
    }

  }

  // broadcast regular correlators to all nodes
  {

    // broadcast number of correlators
    int n = _correlators.size();
    comm.Bcast(&n, 1, MPI::INT, 0);
    if (myrank != 0)
      _correlators.resize(n);

    // broadcast non-ti correlators
    for (typename std::vector< pcps::Correlator<S> >::iterator corr = _correlators.begin(); corr != _correlators.end(); corr++) {
      bool is_ti = corr->ti();
      comm.Bcast(&is_ti, 1, MPI::BOOL, 0);
      if (!is_ti)
        corr->broadcast(comm, myrank, 0);
    }

    // broadcast translationally invariant correlators
    for (typename std::vector< pcps::Correlator<S> >::iterator corr = _correlators.begin(); corr != _correlators.end(); corr++) {
      bool is_ti = corr->ti();
      comm.Bcast(&is_ti, 1, MPI::BOOL, 0);
      if (is_ti) {
        int nsites;
        if (myrank == 0)
          nsites = corr->site_vec().size();
        if (myrank == 0 && nsites == 0)
          throw pcps::Exception("encountered translationally invariant correlator with no sites in pcps::Wavefunction::init_correlators");
        comm.Bcast(&nsites, 1, MPI::INT, 0);
        std::vector<int> sites(nsites);
        if (myrank == 0)
          sites = corr->site_vec();
        comm.Bcast(&sites.at(0), nsites, MPI::INT, 0);
        int parent_index;
        if (myrank == 0)
          parent_index = this->get_ti_parent_index(corr->ti_sites());
        comm.Bcast(&parent_index, 1, MPI::INT, 0);
        if (myrank != 0) {
          pcps::Correlator<S> ti_corr(_correlators.at(parent_index), sites);
          *corr = ti_corr;
        }
      }
    }

  }

  // broadcast perturbative correlators to all nodes
  {
    int n = _pert_corrs.size();
    comm.Bcast(&n, 1, MPI::INT, 0);
    if (myrank != 0)
      _pert_corrs.resize(n);
    for (typename std::vector< pcps::Correlator<S> >::iterator corr = _pert_corrs.begin(); corr != _pert_corrs.end(); corr++)
      corr->broadcast(comm, myrank, 0);
  }

  // compute the positions of each correlator's variables in the ordering of all wavefunction variables
  {

    // assign all offsets to -1 (fixed correlators will keep this value)
    _corr_offsets.assign(_correlators.size(), -1);

    // process non-fixed, non-ti correlators
    int c = 0;
    int offset = 0;
    for (typename std::vector<pcps::Correlator<S> >::const_iterator corr = _correlators.begin(); corr != _correlators.end(); corr++, c++)
      if (!corr->fixed() && !corr->ti()) {
        _corr_offsets.at(c) = offset;
        offset += corr->nelements();
      }

    // process ti correlators (images of translationally invariant correlators)
    c = 0;
    for (typename std::vector<pcps::Correlator<S> >::const_iterator corr = _correlators.begin(); corr != _correlators.end(); corr++, c++)
      if (!corr->fixed() && corr->ti())
        _corr_offsets.at(c) = _corr_offsets.at( this->get_ti_parent_index(corr->ti_sites()) );

  }

}

//-------------------------------------------------------------------------------
// pcps::Wavefunction::init_orb_coeff -- initialize the orbital coefficients
//-------------------------------------------------------------------------------

template <class S>
void pcps::Wavefunction<S>::init_orb_coeff(const pcps::Input & userinp, const pcps::Hamiltonian<S> & ham) {

  // get MPI info
  const MPI::Comm & comm = MPI::COMM_WORLD;
  const int nproc = comm.Get_size();
  const int myrank = comm.Get_rank();

  // size the orbital coefficient array and fill it with zeros
  _restricted_orb_coeffs.assign(userinp.n_rest_orb_coeff(), pcps::zero<S>());

  // get orbital coefficients on root node
  if (myrank == 0) {

    // get the initial orbitals from the binary file
    if (userinp.orb_source() == pcps::keyword::Orb_Source_Binary) {

      // get the name of the binary file
      const std::string filename = userinp.workdir() + "pcps_wfn.bin";

      // open the binary file for reading
      std::ifstream ifs(filename.c_str(), std::ios::in | std::ios::binary);

      // check that the file is open
      if (!ifs.is_open())
        throw pcps::Exception( (boost::format("failed to open file \"%s\" in pcps::Wavefunction::init_orb_coeff") % filename).str() );

      // skip the correlators
      int n;
      ifs.read((char *)&n, sizeof(int));
      for (int i = 0; i < n; i++) {
        pcps::Correlator<S> temp_corr;
        temp_corr.read_binary(ifs);
      }

      // read the number of orbital coefficients
      ifs.read((char *)&n, sizeof(int));

      // check that the number of coefficients is correct
      if (n != _restricted_orb_coeffs.size())
        throw pcps::Exception( (boost::format("wrong number of orbital coefficients in the file \"%s\" ") % filename).str() );

      // read the orbital coefficients
      if (_restricted_orb_coeffs.size() > 0)
        ifs.read((char *)&_restricted_orb_coeffs.at(0), _restricted_orb_coeffs.size() * sizeof(S));

      // close the file
      ifs.close();

    // get the initial orbitals from a text file
    } else if (userinp.orb_source() == pcps::keyword::Orb_Source_Text) {

      // get the name of the text file
      const std::string filename = userinp.workdir() + "pcps_orb_coeff.txt";

      // open the file for reading
      std::ifstream ifs(filename.c_str());

      // check that the file is open
      if (!ifs.is_open())
        throw pcps::Exception( (boost::format("failed to open file \"%s\" in pcps::Wavefunction::init_orb_coeff") % filename).str() );

      // read the orbital coefficients
      for (int i = 0; i < _restricted_orb_coeffs.size(); i++)
        ifs >> _restricted_orb_coeffs.at(i);

      // close the file
      ifs.close();

    // leave the zeros in the orbital coefficients
    } else if (userinp.orb_source() == pcps::keyword::Orb_Source_None) {

    // unknown source for the orbital coefficients
    } else
      throw pcps::Exception( (boost::format("unknown orbital source \"%s\" in pcps::Wavefunction::init_orb_coeff") % userinp.orb_source()).str() );

    // if requested, use a hartree fock calculation to compute the orbitals (using the orbitals from above as the initial guess)
    if (userinp.max_hf_iter() > 0) {

      // full fermions
      if (userinp.particle_type() == pcps::keyword::Part_Fermion)
        pcps::hartree_fock_restricted(userinp, ham, _restricted_orb_coeffs);

      // unsupported particle type
      else
        throw pcps::Exception( (boost::format("cannot run hartree fock calculation for particle type \"%s\"") % userinp.particle_type()).str() );

    }

    // if requested, randomize the orbitals a little
    if (userinp.slightly_random_orbs())
      for (int i = 0; i < _restricted_orb_coeffs.size(); i++) {
        _restricted_orb_coeffs.at(i) += 0.03 * ( pcps::random_number<double>() - 0.5 ) * pcps::unity<S>();
        _restricted_orb_coeffs.at(i) *= 1.0 + 0.05 * ( pcps::random_number<double>() - 0.5 );
      }

  }

  // broadcast the orbital coefficients to all nodes
  if (_restricted_orb_coeffs.size() > 0)
    comm.Bcast(&_restricted_orb_coeffs.at(0), _restricted_orb_coeffs.size() * sizeof(S), MPI::CHAR, 0);

  // initialize 1-body rdm in full alpha beta form (needed for projected optimization and energy)
  if (userinp.optimization() == pcps::keyword::Opt_Projected && userinp.ref_type() == pcps::keyword::Ref_RHF) {

    // spinless fermions
    if (userinp.particle_type() == pcps::keyword::Part_Spinless_Fermion) {

      throw pcps::Exception((boost::format("initialization of hf_rdm not currently implemented for spinless fermions")).str());

    // full fermions
    } else if (userinp.particle_type() == pcps::keyword::Part_Fermion) {

      if (userinp.nparticles() % 2 != 0)
        throw pcps::Exception((boost::format("initialization of hf_rdm for fermions requires the number of particles to be even")).str());

      _hf_rdm.assign(_nsites*_nsites, pcps::zero<S>());

      for (int r = 0; r < _nsites/2; r++)
      for (int s = 0; s < _nsites/2; s++) {
        S x = pcps::zero<S>();
        for (int i = 0; i < userinp.nparticles()/2; i++)
          x += pcps::conj( _restricted_orb_coeffs[i*_nsites/2+r] ) * _restricted_orb_coeffs[i*_nsites/2+s];
        _hf_rdm[ (2*r+0) * _nsites + (2*s+0) ] = x;
        _hf_rdm[ (2*r+1) * _nsites + (2*s+1) ] = x;
      }

    }

  }

}

//-------------------------------------------------------------------------------
// pcps::Wavefunction::init_agp_weights -- initialize the agp weights
//-------------------------------------------------------------------------------

template <class S>
void pcps::Wavefunction<S>::init_agp_weights(const pcps::Input & userinp) {

  // get MPI info
  const MPI::Comm & comm = MPI::COMM_WORLD;
  const int nproc = comm.Get_size();
  const int myrank = comm.Get_rank();

  // size the agp weight array and fill it with zeros
  _agp_weights.assign(userinp.n_agp_weights(), pcps::zero<S>());

  // get weights on root node
  if (myrank == 0) {

    // get the initial weights from the binary file
    if (userinp.agp_source() == pcps::keyword::AGP_Source_Binary) {

      // get the name of the binary file
      const std::string filename = userinp.workdir() + "pcps_wfn.bin";

      // open the binary file for reading
      std::ifstream ifs(filename.c_str(), std::ios::in | std::ios::binary);

      // check that the file is open
      if (!ifs.is_open())
        throw pcps::Exception( (boost::format("failed to open file \"%s\" in pcps::Wavefunction::init_agp_weights") % filename).str() );

      // skip the correlators
      int n;
      ifs.read((char *)&n, sizeof(int));
      for (int i = 0; i < n; i++) {
        pcps::Correlator<S> temp_corr;
        temp_corr.read_binary(ifs);
      }

      // skip the orbital coefficients
      ifs.read((char *)&n, sizeof(int));
      ifs.ignore( n * sizeof(S) );

      // read the number of agp weights
      ifs.read((char *)&n, sizeof(int));

      // check that the number of agp weights is correct
      if (n != _agp_weights.size())
        throw pcps::Exception( (boost::format("wrong number of agp weights in the file \"%s\" ") % filename).str() );

      // read the agp weights
      if (_agp_weights.size() > 0)
        ifs.read((char *)&_agp_weights.at(0), _agp_weights.size() * sizeof(S));

      // close the file
      ifs.close();

    // get the initial weights from a text file
    } else if (userinp.agp_source() == pcps::keyword::AGP_Source_Text) {

      // get the name of the text file
      const std::string filename = userinp.workdir() + "pcps_agp_weights.txt";

      // open the file for reading
      std::ifstream ifs(filename.c_str());

      // check that the file is open
      if (!ifs.is_open())
        throw pcps::Exception( (boost::format("failed to open file \"%s\" in pcps::Wavefunction::init_agp_weights") % filename).str() );

      // read the agp weights
      for (int i = 0; i < _agp_weights.size(); i++)
        ifs >> _agp_weights.at(i);

      // close the file
      ifs.close();

    // use AGP weights that correspond to the restricted Hartree-Fock wavefunction
    } else if (userinp.agp_source() == pcps::keyword::AGP_Source_RHF) {

      // check that the particles are fermions
      if ( userinp.particle_type() != pcps::keyword::Part_Fermion )
        throw pcps::Exception( (boost::format("choosing agp_source as \"%s\" requires particle_type \"%s\"")
                                % pcps::keyword::AGP_Source_RHF
                                % pcps::keyword::Part_Fermion
                               ).str() );

      // set the occupied orbital weights to one and the unoccupied weights to zero
      assert( userinp.nparticles() % 2 == 0 );
      for (int i = 0; i < userinp.nparticles()/2; i++)
        _agp_weights.at(i) = pcps::unity<S>();
      for (int i = userinp.nparticles()/2; i < _agp_weights.size(); i++)
        _agp_weights.at(i) = pcps::zero<S>();

    // leave the zeros in the agp weights
    } else if (userinp.agp_source() == pcps::keyword::AGP_Source_None) {

    // unknown source for the agp weights
    } else
      throw pcps::Exception( (boost::format("unknown agp source \"%s\" in pcps::Wavefunction::init_agp_weights") % userinp.agp_source()).str() );

  }

  // broadcast the agp weights to all nodes
  if (_agp_weights.size() > 0)
    comm.Bcast(&_agp_weights.at(0), _agp_weights.size() * sizeof(S), MPI::CHAR, 0);

}

//-------------------------------------------------------------------------------
// pcps::Wavefunction::init_pairing_matrix -- initialize the pairing matrix
//-------------------------------------------------------------------------------

template <class S>
void pcps::Wavefunction<S>::init_pairing_matrix(const pcps::Input & userinp, const pcps::Hamiltonian<S> & ham) {

  // get MPI info
  const MPI::Comm & comm = MPI::COMM_WORLD;
  const int nproc = comm.Get_size();
  const int myrank = comm.Get_rank();

  // size the pairing matrix array and fill it with zeros
  _pairing_matrix.assign(userinp.n_pairing_variables(), pcps::zero<S>());

  // get pairing matrix on root node
  if (myrank == 0) {

    // get the initial weights from the binary file
    if (userinp.pairing_source() == pcps::keyword::Pairing_Source_Binary) {

      // get the name of the binary file
      const std::string filename = userinp.workdir() + "pcps_wfn.bin";

      // open the binary file for reading
      std::ifstream ifs(filename.c_str(), std::ios::in | std::ios::binary);

      // check that the file is open
      if (!ifs.is_open())
        throw pcps::Exception( (boost::format("failed to open file \"%s\" in pcps::Wavefunction::init_pairing_matrix") % filename).str() );

      // skip the correlators
      int n;
      ifs.read((char *)&n, sizeof(int));
      for (int i = 0; i < n; i++) {
        pcps::Correlator<S> temp_corr;
        temp_corr.read_binary(ifs);
      }

      // skip the orbital coefficients
      ifs.read((char *)&n, sizeof(int));
      ifs.ignore( n * sizeof(S) );

      // skip the agp weights
      ifs.read((char *)&n, sizeof(int));
      ifs.ignore( n * sizeof(S) );

      // read the number of pairing matrix elements
      ifs.read((char *)&n, sizeof(int));

      // check that the number of matrix elements is correct
      if (n != _pairing_matrix.size())
        throw pcps::Exception( (boost::format("wrong number of pairing matrix elements in the file \"%s\" ") % filename).str() );

      // read the pairing matrix elements
      if (_pairing_matrix.size() > 0)
        ifs.read((char *)&_pairing_matrix.at(0), _pairing_matrix.size() * sizeof(S));

      // close the file
      ifs.close();

    // get the initial pairing matrix from a text file
    } else if (userinp.pairing_source() == pcps::keyword::Pairing_Source_Text) {

      // get the name of the text file
      const std::string filename = userinp.workdir() + "pcps_pairing_matrix.txt";

      // open the file for reading
      std::ifstream ifs(filename.c_str());

      // check that the file is open
      if (!ifs.is_open())
        throw pcps::Exception( (boost::format("failed to open file \"%s\" in pcps::Wavefunction::init_pairing_matrix") % filename).str() );

      // read the pairing matrix
      for (int i = 0; i < _pairing_matrix.size(); i++)
        ifs >> _pairing_matrix.at(i);

      // close the file
      ifs.close();

    // create a random pairing matrix
    } else if (userinp.pairing_source() == pcps::keyword::Pairing_Source_Random) {

      // set each pairing matrix element to a random number with a norm that is not too small
      for (int i = 0; i < _pairing_matrix.size(); i++) {
        do {
          _pairing_matrix.at(i) =   pcps::unity<S>() * ( pcps::random_number<double>() - 0.5 )
                                  + pcps::imaginary_unity<S>() * ( pcps::random_number<double>() - 0.5 );
        } while ( std::abs(_pairing_matrix.at(i)) < 1.0e-1 );
      }

      // symmetrize the random pairing matrix
      if (userinp.ref_type() != pcps::keyword::Ref_Pfaffian) {
        assert( _pairing_matrix.size() == userinp.n_agp_weights() * userinp.n_agp_weights() );
        for (int i = 0; i < userinp.n_agp_weights(); i++)
          for (int j = i+1; j < userinp.n_agp_weights(); j++)
            _pairing_matrix.at( userinp.n_agp_weights() * i + j ) = _pairing_matrix.at( userinp.n_agp_weights() * j + i );
      }

    // use a pairing matrix corresponding to the AGP wavefunction
    } else if (userinp.pairing_source() == pcps::keyword::Pairing_Source_AGP) {

      // get the agp matrix dimension
      const int n = userinp.n_agp_weights();
      assert( _agp_weights.size() == n );
      assert( _restricted_orb_coeffs.size() == n*n );

      // zero the pairing matrix
      pcps::xscal(_pairing_matrix.size(), pcps::zero<S>(), &_pairing_matrix.at(0), 1);

      // pairing matrix for pfaffians
      if (userinp.ref_type() == pcps::keyword::Ref_Pfaffian) {

        assert( _pairing_matrix.size() == 2*n*(2*n-1)/2 );
        for (int p = 0; p < n; p++)
        for (int q = 0; q < n; q++)
        for (int i = 0; i < n; i++) {
          const S x = _agp_weights.at(i) * _restricted_orb_coeffs.at(n*i+p) * _restricted_orb_coeffs.at(n*i+q);
          if (p <= q)
            _pairing_matrix.at(pcps::cmpd_pair_index(2*p+0, 2*q+1, 2*n)) += x;
          else
            _pairing_matrix.at(pcps::cmpd_pair_index(2*q+1, 2*p+0, 2*n)) -= x;
        }

      // pairing matrix for other references
      } else {

        assert( _pairing_matrix.size() == n*n );
        for (int p = 0; p < n; p++)
        for (int q = 0; q < n; q++)
        for (int i = 0; i < n; i++)
          _pairing_matrix.at(n*p+q) += _agp_weights.at(i) * _restricted_orb_coeffs.at(n*i+p) * _restricted_orb_coeffs.at(n*i+q);

      }

    // use a pairing matrix derived from a generalized hartree fock guess
    } else if (userinp.pairing_source() == pcps::keyword::Pairing_Source_GHF) {

      if ( userinp.ref_type() != pcps::keyword::Ref_Pfaffian )
        throw pcps::Exception("GHF pairing source is only supported by pfaffian reference functions");

      const int n = userinp.nsites()/2;

      // allocate alpha and beta fock matrices
      std::vector<S> fock_matrix_a(n*n, pcps::zero<S>());
      std::vector<S> fock_matrix_b(n*n, pcps::zero<S>());

      if (userinp.hamiltonian_type() == pcps::keyword::Ham_Hubbard) {

        assert( userinp.hamiltonian_type() == pcps::keyword::Ham_Hubbard );
        assert( userinp.particle_type() == pcps::keyword::Part_Fermion );
        assert( userinp.lattice_type() == pcps::keyword::Lattice_Cartesian );
        assert( userinp.lattice_length_x() == userinp.lattice_length_y() );
        assert( userinp.lattice_length_z() == 1 );

        // build alpha and beta fock matrices
        for (typename std::vector<typename pcps::OpBase<S>::const_ptr>::const_iterator op = ham.ops().begin(); op != ham.ops().end(); op++) {
          if (typeid(**op) == typeid(pcps::OpRL<S>)) {
            assert( (**op).sites().at(0) % 2 == (**op).sites().at(1) % 2 );
            const int site1 = (**op).sites().at(0)/2;
            const int site2 = (**op).sites().at(1)/2;
            if ((**op).sites().at(0) % 2 == 0 )
              fock_matrix_a.at( n * site1 + site2 ) += pcps::unity<S>() * (**op).coeff();
            else
              fock_matrix_b.at( n * site1 + site2 ) += pcps::unity<S>() * (**op).coeff();
          } else if (typeid(**op) == typeid(pcps::OpNN<S>)) {
            const int site = (**op).sites().at(0)/2;
            const double sign = ( site % userinp.lattice_length_x() % 2 == site / userinp.lattice_length_x() % 2 ? 1.0 : -1.0 );
            const double part_per_site = double(userinp.nparticles()) / double(userinp.nsites()/2);
            fock_matrix_a.at(n * site + site) += pcps::unity<S>() * 0.5 * (**op).coeff() * ( part_per_site - sign * userinp.uhf_bias().at(0) );
            fock_matrix_b.at(n * site + site) += pcps::unity<S>() * 0.5 * (**op).coeff() * ( part_per_site + sign * userinp.uhf_bias().at(0) );
          } else {
            throw pcps::Exception( ( boost::format("unexpected operator type \"%s\" in pcps::Wavefunction::init_pairing_matrix")
                                     % pcps::type_lookup(typeid(**op)) ).str() );
          }
        }

        struct Local_Process_Matrix {

          void operator()(const int n, std::vector<S> & fock_matrix, const std::string & name) {

            // print fock matrix
            std::cout << boost::format("printing fock matrix for %s") % name << std::endl;
            std::cout << std::endl;
            for (int i = 0; i < n; i++) {
              for (int j = 0; j < n; j++)
                std::cout << boost::format(" %30s") % pcps::format_number("%10.4f", fock_matrix.at(i*n+j));
              std::cout << std::endl;
            }
            std::cout << std::endl;

            // diagonalize fock matrix
            int info;
            int lwork = 30*n;
            std::vector<double> w(n, 0.0);
            std::vector<S> work(lwork, pcps::zero<S>());
            std::vector<double> rwork(3*n, 0.0);
            std::vector<S> fock_copy(fock_matrix);
            pcps::xsyev('V', 'U', n, &fock_matrix.at(0), n, &w.at(0), &work.at(0), lwork, &rwork.at(0), info);
            pcps::xconj(n*n, &fock_matrix.at(0), 1);

            // check that evecs are in rows
            pcps::xscal(n, pcps::zero<S>(), &work.at(0), 1);
            for (int i = 0; i < n; i++)
            for (int j = 0; j < n; j++)
              work.at(i) += fock_copy.at(n*i+j) * fock_matrix.at(j);
            for (int i = 0; i < n; i++)
              work.at(i) -= w.at(0) * fock_matrix.at(i);
            std::cout << boost::format("error in first eigenpair = %.2e") % std::abs(pcps::xdotc(n, &work.at(0), 1, &work.at(0), 1)) << std::endl;
            std::cout << std::endl;
            assert( std::abs(pcps::xdotc(n, &work.at(0), 1, &work.at(0), 1)) < 1.0e-10 );

            // print evals and evecs
            std::cout << boost::format("printing fock matrix evals and evecs for %s") % name << std::endl;
            std::cout << std::endl;
            std::cout << boost::format("     evals          evecs") << std::endl;
            std::cout << boost::format(" -------------- ------------------------------");
            for (int i = 0; i < n-1; i++)
              std::cout << boost::format("------------------------------");
            std::cout << std::endl;
            for (int i = 0; i < n; i++) {
              std::cout << boost::format(" %14.6f") % w.at(i);
              for (int j = 0; j < n; j++)
                std::cout << boost::format(" %30s") % pcps::format_number("%10.4f", fock_matrix.at(i*n+j));
              std::cout << std::endl;
            }
            std::cout << std::endl;

          }

        };

        Local_Process_Matrix process_matrix;
        process_matrix(n, fock_matrix_a, "alpha");
        process_matrix(n, fock_matrix_b,  "beta");

      // use uhf code
      } else {
        const int nap = userinp.nparticles()/2;
        S energy;
        pcps::hf::uhf(n, nap, nap, ham, userinp.uhf_bias(), userinp.uhf_bias_decay(), 1.0e-9, userinp.max_hf_iter(),
                      fock_matrix_a, fock_matrix_b, energy);
      }

      // randomize orbitals a little if requested
      if (userinp.slightly_random_orbs())
        for (int i = 0; i < n*n; i++) {
          fock_matrix_a.at(i) += 0.03 * ( pcps::random_number<double>() - 0.5 ) * pcps::unity<S>();
          fock_matrix_a.at(i) *= 1.0 + 0.05 * ( pcps::random_number<double>() - 0.5 );
          fock_matrix_b.at(i) += 0.03 * ( pcps::random_number<double>() - 0.5 ) * pcps::unity<S>();
          fock_matrix_b.at(i) *= 1.0 + 0.05 * ( pcps::random_number<double>() - 0.5 );
        }

      // de-segragate alpha and beta
      for (int i = 0; i < n*n; i++) {
        fock_matrix_a.at(i) += userinp.uhf_desegragation() * fock_matrix_a.at(i) / std::abs(fock_matrix_a.at(i));
        fock_matrix_b.at(i) += userinp.uhf_desegragation() * fock_matrix_b.at(i) / std::abs(fock_matrix_b.at(i));
      }

      // get matrix of general hartree fock orbital coefficients B
      assert( userinp.nparticles() % 2 == 0 );
      const int nap = userinp.nparticles()/2;
      std::vector<S> B(userinp.nsites() * userinp.nparticles(), pcps::zero<S>());
      for (int i = 0; i < n; i++)
      for (int j = 0; j < nap; j++) {
        B.at(userinp.nparticles()*(2*i+0)+j) += fock_matrix_a.at(n*j+i) / std::sqrt(2.0);
        B.at(userinp.nparticles()*(2*i+1)+j) += fock_matrix_a.at(n*j+i) / std::sqrt(2.0);
      }
      for (int i = 0; i < n; i++)
      for (int j = 0; j < nap; j++) {
        B.at(userinp.nparticles()*(2*i+0)+j+nap) += fock_matrix_b.at(n*j+i) / std::sqrt(2.0);
        B.at(userinp.nparticles()*(2*i+1)+j+nap) -= fock_matrix_b.at(n*j+i) / std::sqrt(2.0);
      }

      // construct the matrix A, used to convert B into the pairing matrix
      std::vector<S> A(userinp.nparticles() * userinp.nparticles(), pcps::zero<S>());
      for (int i = 0; i < userinp.nparticles()/2; i++) {
        A.at( userinp.nparticles() * (2*i+0) + (2*i+1) ) =  pcps::unity<S>();
        A.at( userinp.nparticles() * (2*i+1) + (2*i+0) ) = -pcps::unity<S>();
      }

      //// construct the matrix A, used to convert B into the pairing matrix
      //std::vector<S> A(userinp.nparticles() * userinp.nparticles(), pcps::unity<S>());
      //for (int i = 0; i < userinp.nparticles(); i++)
      //  A.at(userinp.nparticles()*i+i) = pcps::zero<S>();
      //for (int i = 0; i < userinp.nparticles(); i++)
      //for (int j = 0; j < i; j++)
      //  A.at(userinp.nparticles()*i+j) *= -1.0;

      // construct the pairing matrix
      for (int i = 0; i < userinp.nsites(); i++)
      for (int j = i+1; j < userinp.nsites(); j++) {
        S x = pcps::zero<S>();
        for (int k = 0; k < userinp.nparticles(); k++)
        for (int l = 0; l < userinp.nparticles(); l++)
          x += B.at(userinp.nparticles()*i+k) * A.at(userinp.nparticles()*k+l) * B.at(userinp.nparticles()*j+l);
        _pairing_matrix.at(pcps::cmpd_pair_index(i, j, userinp.nsites())) = x; //pcps::round(x, 6);
      }

    // leave the zeros in the pairing matrix
    } else if (userinp.pairing_source() == pcps::keyword::Pairing_Source_None) {

    // unknown source for the pairing matrix
    } else
      throw pcps::Exception( (boost::format("unknown pairing source \"%s\" in pcps::Wavefunction::init_pairing_matrix") % userinp.pairing_source()).str() );

  }

  // broadcast the pairing matrix to all nodes
  if (_pairing_matrix.size() > 0)
    comm.Bcast(&_pairing_matrix.at(0), _pairing_matrix.size() * sizeof(S), MPI::CHAR, 0);

}

//-----------------------------------------------------------------------------------
//  pcps::Wavefunction::init_touch_vecs -- initializes the vectors of correlators
//                                         telling which correlators touch which sites
//-----------------------------------------------------------------------------------

template <class S> void pcps::Wavefunction<S>::init_touch_vecs() {

  // clear the outer vectors
  _1site_touching_lists.clear();
  _2site_touching_lists.clear();
  _4site_touching_lists.clear();

  // resize the outer vectors
  _1site_touching_lists.resize(_nsites);
  _2site_touching_lists.resize(_nsites * _nsites);
  //_4site_touching_lists.resize(_nsites * _nsites * _nsites * _nsites);

  // allocate the inner vectors
  {

    // single sites
    for (int i = 0; i < _nsites; i++) {
      int n = 0;
      for (int k = 0; k < _correlators.size(); k++)
        if (std::count(_correlators.at(k).site_vec().begin(), _correlators.at(k).site_vec().end(), i) > 0)
          n++;
      _1site_touching_lists.at(i) = boost::shared_ptr<std::vector<pcps::Correlator<S> *> >( new std::vector<pcps::Correlator<S> *>(n) );
    }

    // pairs of sites
    for (int i = 0; i < _nsites; i++)
    for (int j = i; j < _nsites; j++) {
      int n = 0;
      for (int k = 0; k < _correlators.size(); k++)
        if (    std::count(_correlators.at(k).site_vec().begin(), _correlators.at(k).site_vec().end(), i) > 0
             || std::count(_correlators.at(k).site_vec().begin(), _correlators.at(k).site_vec().end(), j) > 0 )
          n++;
      _2site_touching_lists.at(i*_nsites+j) = boost::shared_ptr<std::vector<pcps::Correlator<S> *> >( new std::vector<pcps::Correlator<S> *>(n) );
    }

    //// quadruplets of sites
    //for (int i = 0; i < _nsites; i++)
    //for (int j = i; j < _nsites; j++)
    //for (int p = j; p < _nsites; p++)
    //for (int q = p; q < _nsites; q++) {
    //  int n = 0;
    //  for (int k = 0; k < _correlators.size(); k++)
    //    if ( std::count(_correlators.at(k).site_vec().begin(), _correlators.at(k).site_vec().end(), i) > 0 ||
    //         std::count(_correlators.at(k).site_vec().begin(), _correlators.at(k).site_vec().end(), j) > 0 ||
    //         std::count(_correlators.at(k).site_vec().begin(), _correlators.at(k).site_vec().end(), p) > 0 ||
    //         std::count(_correlators.at(k).site_vec().begin(), _correlators.at(k).site_vec().end(), q) > 0 )
    //      n++;
    //  _4site_touching_lists.at(_nsites*(_nsites*(_nsites*i+j)+p)+q) =
    //    boost::shared_ptr<std::vector<pcps::Correlator<S> *> >( new std::vector<pcps::Correlator<S> *>(n) );
    //}

  }

  // populate the inner vectors
  {

    // single sites
    for (int i = 0; i < _nsites; i++) {
      int n = 0;
      for (int k = 0; k < _correlators.size(); k++)
        if (std::count(_correlators.at(k).site_vec().begin(), _correlators.at(k).site_vec().end(), i) > 0)
          _1site_touching_lists.at(i)->at(n++) = &_correlators.at(k);
    }

    // pairs of sites
    for (int i = 0; i < _nsites; i++)
    for (int j = i; j < _nsites; j++) {
      int n = 0;
      for (int k = 0; k < _correlators.size(); k++)
        if (    std::count(_correlators.at(k).site_vec().begin(), _correlators.at(k).site_vec().end(), i) > 0
             || std::count(_correlators.at(k).site_vec().begin(), _correlators.at(k).site_vec().end(), j) > 0 )
          _2site_touching_lists.at(i*_nsites+j)->at(n++) = &_correlators.at(k);
    }

    //// quadruplets of sites
    //for (int i = 0; i < _nsites; i++)
    //for (int j = i; j < _nsites; j++)
    //for (int p = j; p < _nsites; p++)
    //for (int q = p; q < _nsites; q++) {
    //  int n = 0;
    //  for (int k = 0; k < _correlators.size(); k++)
    //    if ( std::count(_correlators.at(k).site_vec().begin(), _correlators.at(k).site_vec().end(), i) > 0 ||
    //         std::count(_correlators.at(k).site_vec().begin(), _correlators.at(k).site_vec().end(), j) > 0 ||
    //         std::count(_correlators.at(k).site_vec().begin(), _correlators.at(k).site_vec().end(), p) > 0 ||
    //         std::count(_correlators.at(k).site_vec().begin(), _correlators.at(k).site_vec().end(), q) > 0 )
    //      _4site_touching_lists.at(_nsites*(_nsites*(_nsites*i+j)+p)+q)->at(n++) = &_correlators.at(k);
    //}

  }

}

//-----------------------------------------------------------------------------------
//  pcps::Wavefunction::dump_binary -- dumps the wavefunction's correlators to a
//                                     binary file
//-----------------------------------------------------------------------------------

template <class S>
void pcps::Wavefunction<S>::dump_binary(const pcps::Input & userinp, const std::string & fname) const {

  // get the name of the binary file
  const std::string filename = userinp.workdir() + fname;

  // open the binary file for writing
  std::ofstream ofs(filename.c_str(), std::ios::out | std::ios::binary);

  // check that the file is open
  if (!ofs.is_open())
    throw pcps::Exception( (boost::format("failed to open file \"%s\" in pcps::Wavefunction::dump_binary") % filename).str() );

  // write the number of correlators (excluding translationally invariant images)
  int n = 0;
  for (typename std::vector< pcps::Correlator<S> >::const_iterator corr = _correlators.begin(); corr != _correlators.end(); corr++)
    if (!corr->ti())
      n++;
  ofs.write((const char *)&n, sizeof(int));

  // write each correlator
  for (typename std::vector< pcps::Correlator<S> >::const_iterator corr = _correlators.begin(); corr != _correlators.end(); corr++)
    if (!corr->ti())
      corr->dump_binary(ofs);

  // write orbital coefficients
  n = _restricted_orb_coeffs.size();
  ofs.write((const char *)&n, sizeof(int));
  if (n > 0)
    ofs.write((const char *)&_restricted_orb_coeffs.at(0), _restricted_orb_coeffs.size() * sizeof(S));

  // write AGP weights
  n = _agp_weights.size();
  ofs.write((const char *)&n, sizeof(int));
  if (n > 0)
    ofs.write((const char *)&_agp_weights.at(0), _agp_weights.size() * sizeof(S));

  // write pairing matrix
  n = _pairing_matrix.size();
  ofs.write((const char *)&n, sizeof(int));
  if (n > 0)
    ofs.write((const char *)&_pairing_matrix.at(0), _pairing_matrix.size() * sizeof(S));

  // close the file
  ofs.close();

}

// explicitly instantiate the template
template void pcps::Wavefunction< double               >::dump_binary(const pcps::Input &, const std::string &) const;
template void pcps::Wavefunction< std::complex<double> >::dump_binary(const pcps::Input &, const std::string &) const;

//-----------------------------------------------------------------------------------
//  pcps::Wavefunction::to_vector -- convert the correlator variables to a vector
//-----------------------------------------------------------------------------------

template <class S> void pcps::Wavefunction<S>::to_vector(const pcps::Input & userinp, double * const v) const {

  class local_value_setter {
    public:
      static void set_value(const double x, double * const v, const size_t v_ind, const size_t n) {
        v[v_ind] = x;
      }
      static void set_value(const std::complex<double> x, double * const v, const size_t v_ind, const size_t n) {
        v[v_ind]   = x.real();
        v[v_ind+n] = x.imag();
      }
  };

  const size_t n = this->nvars();
  size_t v_ind = 0;
  for (typename std::vector<pcps::Correlator<S> >::const_iterator corr = _correlators.begin(); corr != _correlators.end(); corr++) {
    if (corr->fixed() || corr->ti()) continue;
    for (pcps::Cluster occs(userinp.nsites(), userinp.configs_per_site(), corr->site_vec()); !occs.finished(); occs++, v_ind++)
      local_value_setter::set_value((*corr)[occs.cmpd()], v, v_ind, n);
  }

}

// explicitly instantiate the template
template void pcps::Wavefunction< double               >::to_vector(const pcps::Input & userinp, double * const v) const;
template void pcps::Wavefunction< std::complex<double> >::to_vector(const pcps::Input & userinp, double * const v) const;

//-----------------------------------------------------------------------------------
//  pcps::Wavefunction::from_vector -- get the correlator variables from a vector
//-----------------------------------------------------------------------------------

template <class S> void pcps::Wavefunction<S>::from_vector(const pcps::Input & userinp, const double * const v) {

  class local_value_setter {
    public:
      static void set_value(double & x, const double * const v, const size_t v_ind, const size_t n) {
        x = v[v_ind];
      }
      static void set_value(std::complex<double> & x, const double * const v, const size_t v_ind, const size_t n) {
        x = std::complex<double>(v[v_ind], v[v_ind+n]);
      }
  };

  const size_t n = this->nvars();
  size_t v_ind = 0;
  for (typename std::vector<pcps::Correlator<S> >::iterator corr = _correlators.begin(); corr != _correlators.end(); corr++) {
    if (corr->fixed() || corr->ti()) continue;
    for (pcps::Cluster occs(userinp.nsites(), userinp.configs_per_site(), corr->site_vec()); !occs.finished(); occs++, v_ind++)
      local_value_setter::set_value((*corr)[occs.cmpd()], v, v_ind, n);
  }

}

// explicitly instantiate the template
template void pcps::Wavefunction< double               >::from_vector(const pcps::Input &, const double * const);
template void pcps::Wavefunction< std::complex<double> >::from_vector(const pcps::Input &, const double * const);

//-------------------------------------------------------------------------------
// pcps::Wavefunction::print -- prints the wavefunction
//-------------------------------------------------------------------------------

template <class S> void pcps::Wavefunction<S>::print(const pcps::Input & userinp) {

  // print greeting
  std::cout << boost::format("printing wavefunction") << std::endl;
  std::cout << std::endl;

  // print agp weights and orbital coefficients
  if ( userinp.ref_type() == pcps::keyword::Ref_RHF || userinp.ref_type() == pcps::keyword::Ref_AGP ) {
    std::cout << boost::format("%28s    Orbital Coefficients") % "AGP Weight" << std::endl;
    std::cout << boost::format("----------------------------    ");
    for (int j = 0; j < userinp.n_agp_weights(); j++)
      std::cout << boost::format("------------------------------");
    std::cout << std::endl;
    for (int i = 0; i < userinp.n_agp_weights(); i++) {
      std::cout << boost::format("%28s    ") % pcps::format_number("%10.6f", _agp_weights.at(i));
      for (int j = 0; j < userinp.n_agp_weights(); j++)
        std::cout << boost::format("  %28s") % pcps::format_number("%10.6f", _restricted_orb_coeffs.at(userinp.n_agp_weights() * i + j));
      std::cout << std::endl;
    }
    std::cout << std::endl;
  }

  // print pairing matrix for pairing reference
  if ( userinp.ref_type() == pcps::keyword::Ref_Pairing ) {
    std::cout << boost::format("Pairing Matrix") << std::endl;
    for (int j = 0; j < userinp.n_agp_weights(); j++)
      std::cout << boost::format("------------------------------");
    std::cout << std::endl;
    for (int i = 0; i < userinp.n_agp_weights(); i++) {
      for (int j = 0; j < userinp.n_agp_weights(); j++)
        std::cout << boost::format("  %28s") % pcps::format_number("%10.6f", _pairing_matrix.at(userinp.n_agp_weights() * i + j));
      std::cout << std::endl;
    }
    std::cout << std::endl;
  }

  // print pairing matrix for pfaffian reference
  if ( userinp.ref_type() == pcps::keyword::Ref_Pfaffian ) {
    std::cout << boost::format("Pairing Matrix") << std::endl;
    for (int j = 0; j < userinp.nsites(); j++)
      std::cout << boost::format("------------------------------");
    std::cout << std::endl;
    for (int i = 0; i < userinp.nsites(); i++) {
      for (int j = 0; j < userinp.nsites(); j++) {
        if (i == j)
          std::cout << boost::format("  %28s") % "";
          //std::cout << boost::format("  %28s") % pcps::format_number("%10.6f", pcps::zero<S>());
        else if (i < j)
          std::cout << boost::format("  %28s") % pcps::format_number("%10.6f",  _pairing_matrix.at(pcps::cmpd_pair_index(i,j,userinp.nsites())));
        else
          std::cout << boost::format("  %28s") % "";
          //std::cout << boost::format("  %28s") % pcps::format_number("%10.6f", -_pairing_matrix.at(pcps::cmpd_pair_index(j,i,userinp.nsites())));
      }
      std::cout << std::endl;
    }
    std::cout << std::endl;
  }

  // print correlators
  for (typename std::vector< pcps::Correlator<S> >::const_iterator corr = _correlators.begin(); corr != _correlators.end(); corr++) {
    if (corr->ti()) continue; // do not print duplicates of translationally invariant correlators
    pcps::print_correlator(*corr);
    std::cout << std::endl;
  }

}

// explicitly instantiate the template
template void pcps::Wavefunction< double               >::print(const pcps::Input &);
template void pcps::Wavefunction< std::complex<double> >::print(const pcps::Input &);

//-------------------------------------------------------------------------------
// pcps::Wavefunction::scale_variables -- Function to scale the wavefunction's
//                                        variables so that the largest norm of
//                                        any variable is one.
//-------------------------------------------------------------------------------

template <class S> void pcps::Wavefunction<S>::scale_variables() {
  // scale correlator variables
  for (typename std::vector< pcps::Correlator<S> >::iterator corr = _correlators.begin(); corr != _correlators.end(); corr++)
    if (!corr->fixed() && !corr->ti())
      corr->normalize();
  // scale the AGP weight and orbital coefficient variables
  if (_agp_weights.size() > 0) {
    assert( _agp_weights.size() * _agp_weights.size() == _restricted_orb_coeffs.size() );
    for (int i = 0; i < _agp_weights.size(); i++) {
      S orb_norm = pcps::zero<S>();
      for (int p = 0; p < _agp_weights.size(); p++)
        orb_norm += pcps::square_norm(_restricted_orb_coeffs[i*_agp_weights.size()+p]);
      //std::cout << boost::format("%.4e") % orb_norm << std::endl;
      _agp_weights[i] *= orb_norm;
      pcps::xscal(_agp_weights.size(), pcps::unity<S>() / std::sqrt(orb_norm), &_restricted_orb_coeffs[i*_agp_weights.size()], 1);
    }
    S max_element = pcps::zero<S>();
    for (int i = 0; i < _agp_weights.size(); i++)
      if (pcps::square_norm(max_element) < pcps::square_norm(_agp_weights[i]))
        max_element = _agp_weights[i];
    pcps::xscal(_agp_weights.size(), 1.0 / max_element, &_agp_weights[0], 1);
  }
  // scale the pairing matrix
  if (_pairing_matrix.size() > 0) {
    double m = 0.0;
    for (int i = 0; i < _pairing_matrix.size(); i++)
      m = std::max(m, std::abs(_pairing_matrix[i]));
    pcps::xscal(_pairing_matrix.size(), pcps::unity<S>() / m, &_pairing_matrix[0], 1);
  }
}

// explicitly instantiate the template
template void pcps::Wavefunction< double               >::scale_variables();
template void pcps::Wavefunction< std::complex<double> >::scale_variables();

//-------------------------------------------------------------------------------
// pcps::Wavefunction::Wavefunction -- construct wavefunction from user input
//-------------------------------------------------------------------------------

template <class S>
pcps::Wavefunction<S>::Wavefunction(const pcps::Input & userinp, const pcps::Hamiltonian<S> & ham)
  : _nsites(userinp.nsites())
{

  // get MPI info
  const MPI::Comm & comm = MPI::COMM_WORLD;
  const int nproc = comm.Get_size();
  const int myrank = comm.Get_rank();

  // initialize correlators
  this->init_correlators(userinp);

  // initialize orbital coefficients
  this->init_orb_coeff(userinp, ham);

  // initialize the agp weights
  this->init_agp_weights(userinp);

  // initialize the pairing matrix
  this->init_pairing_matrix(userinp, ham);

  // initialize the lists of which correlators touch which sites
  this->init_touch_vecs();

  // scale the wavefunction's varables to keep norms regular
  this->scale_variables();

}

// explicitly instantiate the template
template pcps::Wavefunction< double               >::Wavefunction(const pcps::Input &, const pcps::Hamiltonian< double               > &);
template pcps::Wavefunction< std::complex<double> >::Wavefunction(const pcps::Input &, const pcps::Hamiltonian< std::complex<double> > &);
